#region Using directives

using System;
using System.Collections;

#endregion

namespace ASharpMachine.SOM
{
	
	/*
	Rather than recalculate the distribution used for the
	changes in the map, we'll do it once and store it.

	*/
	class DistributionMask
	{
		public double[,] mask;
		public int rows;
		public int columns;
		public int centerRow;
		public int centerColumn;

		
		public static double NormalDistribution(double sigma,double mu,double x)
		{
			double normalizer=1.0/(sigma*Math.Sqrt(Math.PI*2));
			double center=(x-mu)*(x-mu);
			double funct=Math.Exp(-center/(2*sigma*sigma));
			return funct*normalizer;
		}
		
		public DistributionMask(int rows,int columns)
		{
			if ((rows%2==0)||(columns%2==0))
			{
				throw new Exception("rows and columns must be an odd number!");
			}
			this.rows=rows;
			this.columns=columns;
			mask=new double[rows,columns];
		}
		
		public void FillMask()
	    {
			centerRow=(rows/2);
			centerColumn=(columns/2);
			double sigma=1.0;
			//Console.WriteLine("{0} {1}",centerRow,centerColumn);
			for (int aRow=0;aRow<rows;aRow++)
			{
				for (int aColumn=0;aColumn<columns;aColumn++ )
				{
					double dr=centerRow-aRow;
					double dc=centerColumn-aColumn;
					double distance=Math.Sqrt(dr*dr+dc*dc);
					//Console.WriteLine("{0} {1} {2} {3} {4}",aColumn,dc,aRow,dr,distance);
					mask[aRow,aColumn]=DistributionMask.NormalDistribution(sigma,0,distance);
				}
			}
	    }

	
	}
	
	
	
	class Map
	{
		private double[,][] vectors;
		private DistributionMask dm;
		private int rows=0;
		private int columns=0;
		private int vectorSize=0;
		
		
		public static void Main()
		{
			Map aMap=new Map(10,10,4);
			
			
		
			
			
			/*Console.WriteLine("Starting speed test");
			for (int i=0;i<1000 ;i++ )
			{
				int[] location=aMap.FindClosestVector(test);
				if (i%10==0)
				{
					Console.Write(".");
				}
			}
			Console.WriteLine("");
			Console.WriteLine("Speed test done");*/
			double[] test=aMap.VectorAt(5,5);
			Console.WriteLine("{0} {1} {2} {3}",test[0],test[1],test[2],test[3]);
			double[] testVector={3,3,3,3};
			aMap.transform(5,5,testVector,0.4);	
			test=aMap.VectorAt(7,7);
			Console.WriteLine("{0} {1} {2} {3}",test[0],test[1],test[2],test[3]);
			test=aMap.VectorAt(6,6);
			Console.WriteLine("{0} {1} {2} {3}",test[0],test[1],test[2],test[3]);
			test=aMap.VectorAt(4,4);
			Console.WriteLine("{0} {1} {2} {3}",test[0],test[1],test[2],test[3]);


			
			//ShowMatrix(dm.mask,5,5);
		}

		public static void ShowMatrix(double[,] aMatrix,int rows,int columns)
		{
			for (int aRow=0;aRow<rows;aRow++)
			{
				for (int aColumn=0;aColumn<columns ;aColumn++ )
				{
					Console.Write("{0} ",aMatrix[aRow,aColumn]);
				}
				Console.WriteLine();
			}
		
		}
		
		public static double VectorDistance(double[] v1,double[] v2)
		{
			if(v1.Length!=v2.Length)
			{
				throw new Exception("Vectors must be of the same length!");
			}
			else
			{
				double sum=0;
				//do the euclidian formula thing
				for (int i=0;i<v1.Length ;i++ )
				{
					double diff=v1[i]-v2[i];
					sum+=diff*diff;
				}
				return Math.Sqrt(sum);
			}
		}
		
		public static double[] VectorSub(double[] v1, double[] v2)
		{
			if(v1.Length!=v2.Length)
			{
				throw new Exception("Vectors must be of the same length!");
			}
			else
			{
				double[] rv=new double[v1.Length];
				for (int i=0;i<v1.Length ;i++ )
				{
					rv[i]=v1[i]-v2[i];	
				}
				return rv;
			}
			
		}
		
		public static double[] VectorAdd(double[] v1, double[] v2)
		{
			if(v1.Length!=v2.Length)
			{
				throw new Exception("Vectors must be of the same length!");
			}
			else
			{
				double[] rv=new double[v1.Length];
				for (int i=0;i<v1.Length ;i++ )
				{
					rv[i]=v1[i]+v2[i];	
				}
				return rv;
			}
		}
		
		public static double[] VectorMul(double[] v1, double m)
		{
			double[] rv=new double[v1.Length];
			for (int i=0;i<v1.Length ;i++ )
			{
				rv[i]=v1[i]*m;	
			}
			return rv;
		}

		
		public  Map(int rows, int columns, int vectorSize)
		{
			vectors=new double[rows,columns][];
			this.rows=rows;
			this.columns=columns;
			this.vectorSize=vectorSize;
			
			for (int aRow=0;aRow<rows;aRow++)
			{
				for (int aColumn=0;aColumn<columns ;aColumn++ )
				{
					//we are relying on the compiler to put everything
					//into a known state.
					vectors[aRow,aColumn]=new double[vectorSize];
				}
			}
			dm=new DistributionMask(5,5);
			dm.FillMask();
		}
	
		public double[] VectorAt(int row, int column)
		{	
			return vectors[row,column];
		}

		public void SetVectorAt(int row, int column, double[] aVector)
		{
			if (aVector.Length==vectorSize)
			{
				vectors[row,column]=aVector;
			}
			else
			{
				throw new Exception("Presented vector has different length than map!");
			}
		
		}

		public int[] FindClosestVector(double[] aVector)
		{
			/*set the intial distance as the distance from the 
			 (0,0) vector */
			double initialDistance=Map.VectorDistance(aVector,VectorAt(0,0));
			int[] vectorLocation={0,0};
			
			for (int aRow=0;aRow<rows;aRow++)
			{
				for (int aColumn=0;aColumn<columns ;aColumn++ )
				{
					double checkDistance=Map.VectorDistance(aVector,VectorAt(aRow,aColumn));
					if (checkDistance<initialDistance)
					{
						initialDistance=checkDistance;
						vectorLocation[0]=aRow;
						vectorLocation[1]=aColumn;
					}
				}
			}
			return vectorLocation;

		}
	
		/*
		Given aRow and aColumn this method moves the internal vector at that point
		closer to the seed vector by amount rate and the internal distribution mask
		It then moves all of the neighbors closer
		*/
		public void transform(int aRow,int aColumn,double[] seedVector, double rate)
		{
			//The mask is centered on aRow, aColumn, so we need to start from either the
			//edge of the internal map or the edge of the mask, whichever is lesser,
			//then move to the edge of the map or the edge of the mask, whichever is lesser.
			//some fun could be had here if we postulated a torus for the map. (perhaps a delegate)
			int minCol=Math.Max(0,aColumn-dm.centerColumn);
			int maxCol=Math.Min(columns-1,aColumn+dm.centerColumn);
			int minRow=Math.Max(0,aRow-dm.centerRow);
			int maxRow=Math.Min(rows-1,aRow+dm.centerRow);
			int rowOffset=dm.centerRow-aRow;
			int columnOffset=dm.centerColumn-aColumn;
			//Console.WriteLine("{0}->{1} {2}->{3}",minCol,maxCol,minRow,maxRow);
			for (int i=minRow;i<=maxRow;i++)
			{
				for (int j=minCol;j<=maxCol;j++)
				{
					int idm=i+rowOffset;
					int jdm=j+columnOffset;
					//Console.WriteLine("{0} {1}",idm,jdm);
					//need a vector datatype. badly. the following is crufty!
					double[] vectorDelta=Map.VectorMul(Map.VectorMul(Map.VectorSub(seedVector,vectors[i,j]),dm.mask[idm,jdm]),rate);
					vectors[i,j]=Map.VectorAdd(vectors[i,j],vectorDelta);
				}
			}
		}


	}
	
}

